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rate of the Weibel instability for a relativistic Maxwellian distribution and compare 
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dispersion relation for the relativistic tearing mode, making no assumption about the 
thickness of the current sheet, and check the numerical method against the analytical 
expression. 
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1. Introduction 

Dissipation of the energy carried by relativistic plasma outflows is important for the 
physics of pulsar winds and gamma-ray bursts (for recent reviews see [1] and [2]). In 
these objects, the plasma is probably composed of electrons, positrons and protons. 
As well as being in relativistic bulk motion with respect to the observer, the random 
thermal energy of the plasma may also be relativistic, i.e., comparable to the rest mass 
energy of the constituent particles. In this paper we concentrate on two instabilities: 
these are the two-stream or Weibel instability |3] and the tearing modes in a relativistic 
neutral pair plasma current sheet, thought to play a role in the formation process of 
relativistic shocks [3] and the dissipation of magnetic energy in pulsar winds [5]. They 
are investigated by generalising and extending to the inhomogeneous case the method 
presented in [6] . Motivated primarily by the need for code veriflcation, we have derived 
some analytical expressions for the linear growth rates of these instabilities. 

The Weibel instability is very important in astrophysical processes because it is able 
to generate a magnetic fleld by extracting the free energy from an anisotropic momentum 
distribution in an unmagnetised plasma or from the kinetic drift energy. There is an 
extensive literature on the Weibel instability: general conditions for the existence of the 
relativistic Weibel instability for arbitrary distribution functions are discussed in [7] , and 
wave propagation in counter-streaming magnetised nonrelativistic Maxwellian plasmas 
are studied in [H [9] . Dispersion curves have been found in some special cases such as, for 
example, the fully relativistic bi-Maxwellian distribution function, (Yoon [12]), and the 
water-bag distribution function, in which case closed-form analytical expressions can be 
derived not only for the Weibel instability (Yoon [lOJ), but also for the cyclotron maser 
and whistler instabilities (Yoon [llj). However, flnding an analytical expression for the 
dispersion relation for a given equilibrium distribution function is a complicated or even 
impossible task. It involves a four-dimensional integration (3D in momentum space 
and ID in time) of the equilibrium distribution function which is difficult to perform 
in closed form. For this reason, the water-bag distribution is the preferred profile to 
analyse magnetic field generation in fast ignitor scenarios, (Silva et al. [13]) and in 
relativistic shocks, (Wiersma and Achterberg [H], Lyubarsky and Eichler [15]). The 
Weibel instability in a magnetised electron-positron pair plasma has been investigated 
by Yang et al [16] using two model distributions: the water bag, and one with a power- 
law dependence at high energy. A general covariant description has been formulated 
by Melrose [17] and by Schlickeiser [18]. In the present work, we focus on equilibrium 
configurations given by a relativistic Maxwellian distribution function, which allows 
one to reduce the four- dimensional integral to a simple one-dimensional integral, as we 
demonstrate in Section 12.21 The growth rates are then found by solving this equation 
using a single numerical quadrature, and are compared to the results found using our 
extended algorithm in Section 14.11 

In the inhomogeneous case, the stability properties of a nonrelativistic Harris 
current sheet also have a substantial literature, with notable recent studies by 
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Daughton [T^] and Silin et al. [20] • In the relativistic case, the tearing mode instabihty 
has been investigated by Zelenyi & Krasnoselskikh [21], by integrating first order 
perturbations of the relativistic Maxwelhan distribution function along approximate, 
straight-line particle trajectories, in the thick layer limit (in which the Larmor radius 
is much smaller than the thickness of the current sheet). In Section [2.41 we lift these 
restrictions to present new results for the tearing mode instability in a neutral current 
sheet of arbitrary temperature and thickness, and compare these with the results 
found using the generalised algorithm in Section I4.2[ In this work, no assumption 
is made about the thickness of the current sheet, and the particle trajectories are found 
numerically in the background magnetic field. 

The numerical method, which is an extension of our previous algorithm, [6], that 
computes the linear dispersion relation of waves within a Vlasov-Maxwell description, 
is described in Section [3l It is based on the approach of Daughton pjj for non 
relativistic Maxwellians, and involves explicit time integration of particle orbits along 
the unperturbed trajectories. We modify and extend our former code to include 
inhomogeneities in the plasma equilibrium configuration. Moreover, we generalise it 
to a fully relativistic approach, i.e., one that allows for relativistic temperatures as well 
as relativistic drift speeds. 



2. Analytical treatment of the instabilities 

2.1. The Vlasov-Maxwell equations 

For convenience, we first recall the full set of Vlasov-Maxwell equations governing the 
non-linear time evolution of the pair plasma. 

We introduce the standard electromagnetic scalar and vector potentials {(j),A), 
related to the electromagnetic field {E, B) by : 

E-~V,-- (1) 

B = V ^A (2) 
We employ the Lorenz gauge condition by imposing 

<1\Y A + Eono— = Q (3) 

with £0 fJ'OC^ = ^ and c the speed of hght. The relation between potentials and sources 
then reads : 

A^-^^^+£ =0 (4a) 

^ — * 

The source terms represented by the charge p and current j densities, are expressed in 
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terms of the distribution functions of each species, f., by : 



j(f, t) = ^ yyy /,(f, p, t) (s^) 



where 7 = \/l + p'^/m^ (? is the Lorentz factor of a particle. We adopted the usual 
notations, namely {t,r,v,p,ms,qs) for respectively the time, position, 3- velocity, 3- 
momentum, mass and charge of a particle of species s. The time evolution of the 
distribution functions fs is governed by a relativistic Vlasov equation for each species: 

f + + ,.(B%,7a5).|^0 (6) 

The self-consistent non-linear evolution of the plasma is entirely determined by the set 
of equations ([I])-©. 

The linear stability properties of the neutral current sheet are investigated by 
linearising the set of equations ([I])-©. The procedure has been described in [6j, Sec. 3. 
These results are now used to derive linear dispersion relations. 

2.2. Weibel instability 

In the Weibel regime, the electrostatic potential can be neglected. Moreover, only 
the Ay component of the vector potential A comes into play in the eigenvalue problem. 
Note also that the Weibel instability corresponds to low-frequency modes propagating 
along the 2;-axis such that 

I |a;| I <^ c (7) 

with ky = 0, k = (0, ky, k^) being the wavenumber and ui the corresponding 
eigenfrequency. 

Specialising the general derivation of the perturbed charge density and current 
density presented in Sect. 3 of [6J to the Weibel instability, setting (p = A^ = A^ = ky = 
0, the dispersion relation reads 

X / Pyexp[z (fczPzMs - 7^^) ^Vt'c?V= (8) 

The distribution function at equilibrium for each species "s", denoted by fos{r,p), is 
assumed to be a relativistic Maxwellian with constant drift speed ±[7^. The stationary 
distribution function for each species reads: 

Ng is the (constant) particle number density of the plasma, Upg = Ng ql/ma Eq the plasma 
frequency, E = mgC^ ■>/! + p"^ /ml the total energy of a particle, Py the y-component 



of its momentum, Tg = 1 / a/1 — Ps iPs = Ug/c) the Lorentz factor associated with the 
drift motion and K2 the modified Bessel function of the second kind and of order 2. 
The temperature of the gas T, is conveniently normahsed to the rest mass energy of the 
leptons such that 

e. = ^ (10) 

and kB is the Boltzmann constant. The triple integral along the momentum vector can 
be performed analytically in the following way. 

Because no external background electromagnetic field exists at equilibrium, the 
particle trajectories are straight lines. The integration of the equations of motion leads 
to 

p' = p = est (11) 

f' = f+^r' (12) 
rus 

Therefore time and momentum integration can be inverted. The dispersion relation 
thus reduces to 

exp(-r, (E - f/,py)/esmsc2) / n 

The integration along the momentum p can be done analytically with help on the 
following formula, see for instance [22j 



r/ . -\ 1 /"/"/" exp(-A7 ±ia-p) n nKi{J + m? 6 ^ j ,^ 

I{A,a) = — ^ ^d'p = mlc'^^f===i=^{U) 



where the Lorentz factor of a particle is 7 = a/1 + p'^ /ml (? and K\ is the modified 
Bessel function of order 1. By differentiating twice with respect to the y component, 
Oy, of the vector a, we get 

9^/,, 1 /'/'/' Pyexp(-A7±ia-p) 



/p.(A«)- -^(A«) = ^ yyy ^ -'^'v= (is) 



m\ e' 



^ A? + — 4 m\ (? al , 



ml (? Oy 



(A^ + a^)^ 



{A^ + ml a2)3/2 

Applying these formulae to our problem, it is convenient to introduce the following 
quantities 

A{uJ,T) = ^+iuJT (16) 
^ s 

a(r) = - z — Cy H (17) 

Q.m.c m. 
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The function Ipy depends now on u and r via A{u,t) and a(r), Jpy(A(co', r), a(r)), 
assuming that the equihbrium quantities such as species temperatures 9^ and drift 
speeds are prescribed and therefore constant. In the remainder of this section, we will 
denote it for simplicity by /py(ci;,r). The dispersion relation for the Weibel instability, 
Eq. (|T3ll therefore reads 

In the non-relativistic limit, <^ 1, we recover the dispersion relation Eq.(33) of [6]. 
2.3. Relativistic neutral current sheet 

We now turn to the study of the unstable tearing mode in a relativistic neutral current 
sheet of thickness L and asymptotic magnetic field intensity (far from the current 
sheet). The plasma is one-dimensional in the sense that it has only spatial variation in 
the X direction. It consists of counter-streaming electrons and positrons with relativistic 
temperatures Tg evolving in a static external magnetic field aligned with the z-a.xis such 
that [23] 

B,{f) = Bo tanh (|) (19) 

We use Cartesian coordinates, denoted by r = {x,y,z), and the corresponding 
basis (ex,ey, Cz). In equilibrium, there is no electric field, Eq = and the charges 
drift in the y direction at a relativistic velocity Ug. The particle number density for 
each species is 

2 / X 



rio.(r) = A^ssech" (^-j (20) 

The distribution function at equilibrium for each species "s", denoted by fos{r,p), is 
assumed to be a relativistic Maxwellian with constant drift speed if/^. The stationary 
distribution function for each species reads: 

^'^^^^^P^ = I rT?r^^7T7?TT"^P[-r^ - Ugp,)/Qgmgc'] (21) 

47rmf c'^9sA2(l/9s) 

nos(r) is the particle number density in the current sheet, Eq. fl2U]) . and all other 
quantities are the same as those defined for the distribution function in Eq. (Q. The 
stationary Vlasov-Maxwell equations are satisfied provided that 

ANgQsmgC^ = ^ (22) 

Be nisC^ , ^ 

TsUs = -2 W 23 

Eq. fl22l) states the balance between gaseous pressure and magnetic pressure at the centre 
of the current sheet. It is also useful to note that the non-relativistic plasma frequency 
(at the centre of the current sheet) and cyclotron frequency defined, respectively, by 

= ^ (24) 
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= ^ (25) 

are related by 

4c2;J^e, = fiL (26) 

These frequencies are constant (independent of x), and are denoted by a "in order to 
distinguish them from quantities that depend on x. 

2.4- The relativistic tearing mode 

A method similar to that used in Section 12.21 can also be applied to the relativistic 
neutral current sheet. Using the equilibrium distribution function presented in Eq. fl2Tl) 
for the relativistic Harris current sheet, the eigenvalue equation to solve reads 

Kix) - (^kl - A,ix) + -| sech^ (I) A,{x) + 

where prime " denotes second derivative with respect to x. We now follow [23] and 
introduce the variable t = tanh(a:/L), which transforms Eq. fl27j) into the Legendre 
equation. From this, one sees that the dispersion relation is satisfied by purely growing 
modes, Re(ci;) = 0, the tearing modes, whose growth rate is the solution of 



I Ipy{uj,r)dT={k,L + 2){k,L-l) (28) 

J — oo 



c^eiK2{i/es) mlc> J_, 

where we used the "low- frequency" approximation, Eq. ([7]), to neglect the contribution 
of the displacement current. Expression (1281) contains a single one-dimensional integral. 
Finding its solutions is computationally much faster than dealing with the general 
expression in 4 dimensions and does not involve any approximation concerning the 
current sheet thickness. 

In the low-temperature limit, ^ 1 (non-relativistic case), the dispersion relation 
Eq. (1271) reduces to its classical expression given by 



where the plasma dispersion function Z is defined for Im((^) > by 



1 r'^ e" 

Z{0 = ^ / f-dt (30) 



and analytically continued for Im((^) < 0, see for instance Delcroix and Bers [25] 
3. The algorithm 



We first recall the general linear eigenvalue system to be solved, as presented in [6j and 
then present the extended algorithm for inhomogeneous and magnetised plasmas. 
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3.1. The eigenvalue system 

The eigenvalue system is found by solving the equations for the electromagnetic potential 
determined according to the source distribution given by ( l5al) and (ISfej) . Inserting the 
latter expressions into (l4al) and (l46l) . the eigenvalue system reads : 

<P"{x)~[k^-'^^ <P{x) + '^ =0 (31a) 

i"(x) -{k^-^^^ A{x) + /io /(x) = (316) 

For homogeneous plasmas, the system reduces to a 4x4 matrix, which is easily solved, see 
Pj. In the inhomogeneous case, the perturbations should vanish asymptotically, when 
X = ±oo. To solve this inhomogeneous eigenvalue problem, it is convenient to expand 
the unknown quantities on a set of basis functions, each function individually satisfying 
the required boundary conditions. This method is known as a spectral Galerkin method, 
similar to Fourier decomposition for periodic boundary conditions. We will refer to it 
as a spectral decomposition as described below. 

3.2. Spectral decomposition 

The perturbations are projected on an orthonormal set of basis functions, Because 
of the required boundary conditions, namely 0(a; = ±oo) = 0, A{x = ±oo) = 0, the 
Hermite functions are a convenient set onto which to expand the perturbations, |19j . 
They are given by 

M^) = -4==p exp[-xV2] (32) 
a/2" n\ v^vr 

where if„ are the Hermite polynomials, [26]. We introduce the projection operator of a 
function g on a. basis function JF„ defined as 



<g\J^n>= / g{x)Tnix)dx (33) 

J — oo 

The basis functions themself satisfy the orthonormality relation 

< J^nlJ^k >= Snk (34) 

where Snk is the kronecker symbol. The matrix representation of the second derivative, 
useful for projection of Eq. (131 al) and (13161) . is then represented by 

/ + 00 
J^n{x) J^kix) dx 
-oo 

V(n + l)(n + 2) ^ 2n+l ^ , Vn{n-1) ^ 

— On,k+2 Onk H On,k-2 (,<J0J 

For convenience, we introduce the unknown four dimensional vector \[' = {(f),A). Using 
N terms in the expansion of \E', it is written as 

N-l 

^{x) = Y,^kMx) (36) 

fc=0 
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Inserting this expression in Eq. fl31ap and (131 6p and projecting them on each basis 
function JF„, the system reduces to a matrix equation of dimension AN x AN for the 
unknown vectors 

M{uj, k)-^ = (37) 

More exphcitly, let Ck be the AN unknowns. We choose to order the unknowns in the 
following way, = {Ck,Ck+N,Ck+2N,Ck+3N)- The expansion Eq. fl36l) is done such 
that 

7V-1 

= ^Cfc^fc(x) (38 a) 

fe=0 

7V-1 



A,{x) = Y.Ck+NJ'kix) (386) 

A:=0 
Af-1 

Ay{x) = '^Ck+2N J^k{x) (38c) 



A:=0 
Af-1 



Mx) = 5^Cfe+3jv-^fc(x) (38cO 

k=0 

Eq. (1371) is a non-linear eigenvalue problem for the matrix M with eigenvector 
and eigenvalue u. The method of solution has been described in |6]. However, the 
simultaneous search of the eigenvalues and the eigenvectors is very time-consuming 
because the root finding takes place in a 8A^ + 2 dimension parameter space {AN + 1 
complex numbers to find, the eigenvalue u and the Ck)- Therefore, we decide only to 
look for the vanishing of the matrix determinant, detM(co',A;) = 0. In this paper, we 
only present the capabilities of our algorithm and do not discuss in detail the physics 
of the unstable perturbations. For this purpose, it is sufficient to look only for the 
vanishing of the determinant of the matrix M. Nevertheless, we present an example 
of an eigenfunction in Section 14.21 where we examine the convergence properties of the 
summations in Eqs. fl38fflH38dj) - 



4. Results 



4-1- Weibel instability 

As a first check, we compute the dispersion relation for the Weibel instability for 
arbitrary plasma temperature G^. The fully relativistic algorithm is checked against the 
numerical solution to the exact analytical dispersion relation, Eq. ( ITSl) . The evolution of 
the growth rates with increasing plasma temperature at a given drift speed is shown in 
fig[Tl In our previous work, [6], the non-relativistic case was solved with a non-relativistic 
code whereas in the present work, our algorithm deals with arbitrary temperatures 
(especially in the limit ^ IjtM- In particular, the results of the non-relativistic 



I However, this limit needs special care because of the exponentially decreasing of the Bessel 
function K2 in the expressions Eq. (|21|) . implying numerical under- and overflow problems. 
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Weibel instability are recovered with the fully relativistic algorithm. Inspecting fig[Il the 
solutions of our algorithm (denoted by symbols) are in perfect agreement with the growth 
rates given by numerically solving Eq. ( ITSl) (solid lines), for arbitrary temperature. Due 
to the spread in momentum present at finite temperature, the instability is suppressed 
for large wavenumbers k^,. More precisely, when the temperature increases, the largest 
unstable wavenumber decreases. For a given drift speed, the growth rates decrease with 
increasing temperature. 



Weibel instability 




-4-2 2 



Figure 1. Growth rates of the Weibel instabihty for different temperatures 8^, see 
legend, and a given drift speed (3s — O.f . Points are results from the algorithm and 
the solid lines are from the numerical solution to the dispersion relation, Eq. (jlSp . 



We also investigate the effect of the drift speed on the growth rates for a fixed 
temperature. Results are shown in the non-relativistic limit 0^ = 10~^, 10~^, the mildly 
relativistic case = 1, and the ultra-relativistic case 0^ = 10^, for increasing drift 
speed Ps = 0.1/0.3/0.5, respectively red diamonds, green stars and blue squares, fig. [21 

The characteristic cut-off wavenumber, fccut, can be estimated by setting u; = in 
the dispersion relation Eq. (ITSl) . Doing so we find 




^ps V 

Moreover, in the ultrarelativistic limit, for high temperature Gs ^ 1 but small drift 
speeds /5s <^ 1, the growth rates in the small wavenumber limit are simply given by 

7rel = -/9'fczC (40) 
TT 

This result is in agreement with the curves shown in fig. [TJ Note that these growth 
rates become independent of the temperature, as in the non-relativistic limit (low 
temperature), where 

7clas = Pshc (41) 

Therefore, in the ultrarelativistic temperature limit, the growth rates are reduced by a 
factor 4 ^g./vr ^ 1.27/3,. 
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Weibel instability Weibel instability 




Figure 2. Growth rates of the Weibel instabihty for different drift speeds f3s = 
0.1/0.3/0.5, respectively red diamonds, green stars and blue squares, and a given 
temperature, 9s = lO"'^ in fig a), 6^ = 10"^ in fig b), 6s = 1 in fig c), 0s = 10^ in 
fig d). Points are from our algorithm whereas the solid lines are from the numerical 
solution to the dispersion relation, Eq. ([T5)). 



4-2. Tearing mode 

Next, we check our algorithm against the tearing mode instabihties in the relativistic 
Harris current sheet. The eigenvalue problem for the non-relativistic Harris current 
sheet is given by Galeev, chapter 6.2, page 305, in The dispersion relation for the 
non-relativistic tearing modes {ky = 0) is given in p^[28] . The relativistic generalisation 
is given by Eq. ( !28l) . In figure [HI we compare our numerical results with the analytic 
expressions found using a simple root finding algorithm from the dispersion relation for 
tearing modes in order to check the correctness of our algorithm implementation. The 
growth rates are normalised to the non-relativistic cyclotron frequency CIb = gsBo/nis- 
The numerical results in fig [3] are in good agreement with the approximated analytical 
expression, Eq. ( 128|) . In the spectral decomposition method, it is necessary to choose 
an appropriate number of terms in the expansion (136|1 in order to achieve a given 
accuracy. To do this, we check the convergence properties of our method. In Fig. H] 
we show the eigenvalue as a function of A^, starting with = 3 terms and increasing 
N until the eigenvalue reaches a satisfactory precision (in this case 1 % accuracy). The 
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Relativistic Harris sheet 
0.2| 




0.2 0.4 0.6 0.8 1 

k, L 



Figure 3. Growth rates of the tearing modes in the Harris current sheet for p = L. 
Points are from our numerical algorithm and the solid line represents the analytical 
approximation in Eq. (I28p . 



corresponding convergence of the eigenfunction is shown in fig [5l Note that the results 
of fig [3] are shown for = 11. 
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Convergence properties 
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Figure 4. Example of convergence of an eigenvalue when the number of terms 
in the expansion increases. This specific example corresponds to L = 0.85 and 
Logio ©. = -4 of fig S 



5. Conclusion 

We present a generalisation and extension of our previous algorithm |6] to solve 
the linear dispersion relation for relativistic multi-component inhomogeneous and 
magnetised plasmas. The code is validated by comparing the results with two standard 
configurations: the relativistic Weibel instability in a homogeneous plasma, and the 
tearing mode instability in a relativistic neutral Harris sheet. To effect the comparison, 
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Eigenf unction convergence 




X 



Figure 5. Example of convergence of the eigenfunction associated to tlie eigenvalue in 
figSl Only the coefficients associated with the even Hermite functions are different from 
zero as expected from symmetry considerations. The narrowest function corresponds 
to = 3 whereas the widest to N — 15. 

we derived useful analytical expressions, Eq. f|T8l) and Eq. fl28|) . for the dispersion 
relations in these configurations and solved them numerically. We conclude that this 
code is a suitable tool for the study of stability properties of more general configurations 
of interest in gamma-ray burst and pulsar wind theories. 
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